home *** CD-ROM | disk | FTP | other *** search
/ Nebula 2 / Nebula Two.iso / Apps / Astro / astrolog / Source / formulas.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-12  |  30.8 KB  |  1,061 lines
  1. /*
  2. ** Astrolog (Version 4.10) File: formulas.c
  3. **
  4. ** IMPORTANT NOTICE: the graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1994 by Walter D. Pullen
  6. ** (cruiser1@stein.u.washington.edu). Permission is granted to freely
  7. ** use and distribute these routines provided one doesn't sell,
  8. ** restrict, or profit from them in any way. Modification is allowed
  9. ** provided these notices remain with any altered or edited versions of
  10. ** the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 3/19/1994.
  36. */
  37.  
  38. #include "astrolog.h"
  39.  
  40. real MC, Asc, Vtx, RA, OB, X, Y, A;
  41.  
  42. real planet1[TOTAL+1], planet2[TOTAL+1],
  43.   planetalt1[TOTAL+1], planetalt2[TOTAL+1],
  44.   house1[SIGNS+1], house2[SIGNS+1], ret1[TOTAL+1], ret2[TOTAL+1];
  45.  
  46. real FAR *datapointer;
  47.  
  48.  
  49. /*
  50. ******************************************************************************
  51. ** Specific Calculations.
  52. ******************************************************************************
  53. */
  54.  
  55. #ifdef MATRIX
  56. /* Given a month, day, and year, convert it into a single Julian day value, */
  57. /* i.e. the number of days passed since a fixed reference date.             */
  58.  
  59. long MdyToJulian(mon, day, yea)
  60. int mon, day, yea;
  61. {
  62. #ifndef PLACALC
  63.   long im, j;
  64.  
  65.   im = 12*((long)yea+4800)+(long)mon-3;
  66.   j = (2*(im%12) + 7 + 365*im)/12;
  67.   j += (long)day + im/48 - 32083;
  68.   if (j > 2299171)                   /* Take care of dates in */
  69.     j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  70.   return j;
  71. #else
  72.   int greg = TRUE;
  73.  
  74.   if (yea < G2JYEA || (yea == G2JYEA &&
  75.     (mon < G2JMON || (mon == G2JMON && day < 15))))
  76.     greg = FALSE;
  77.   return (long)floor(julday(mon, day, yea, 12.0, greg)+ROUND);
  78. #endif
  79. }
  80.  
  81.  
  82. /* Take a Julian day value, and convert it back into the corresponding */
  83. /* month, day, and year.                                               */
  84.  
  85. void JulianToMdy(JD, mon, day, yea)
  86. real JD;
  87. int *mon, *day, *yea;
  88. {
  89. #ifndef PLACALC
  90.   long L, N, IT, JT, K, IK;
  91.  
  92.   L  = (long)floor(JD+ROUND)+68569L;
  93.   N  = Dvd(4L*L, 146097L);
  94.   L  -= Dvd(146097L*N + 3L, 4L);
  95.   IT = Dvd(4000L*(L+1L), 1461001L);
  96.   L  -= Dvd(1461L*IT, 4L) - 31L;
  97.   JT = Dvd(80L*L, 2447L);
  98.   K  = L-Dvd(2447L*JT, 80L);
  99.   L  = Dvd(JT, 11L);
  100.   JT += 2L - 12L*L;
  101.   IK = 100L*(N-49L) + IT + L;
  102.   *mon = (int)JT; *day = (int)K; *yea = (int)IK;
  103. #else
  104.   int greg = TRUE;
  105.   double tim;
  106.  
  107.   if (JD < 2299171.0)    /* October 15, 1582 */
  108.     greg = FALSE;
  109.   revjul(JD, greg, mon, day, yea, &tim);
  110. #endif
  111. }
  112.  
  113.  
  114. /* This is a subprocedure of CastChart(). Once we have the chart parameters, */
  115. /* calculate a few important things related to the date, i.e. the Greenwich  */
  116. /* time, the Julian day and fractional part of the day, the offset to the    */
  117. /* sidereal, and a couple of other things.                                   */
  118.  
  119. real ProcessInput(var)
  120. int var;
  121. {
  122.   real Off, Ln;
  123.  
  124.   TT = Sgn(TT)*floor(dabs(TT))+FRACT(dabs(TT))*100.0/60.0+DecToDeg(ZZ);
  125.   OO = DecToDeg(OO);
  126.   AA = MIN(AA, 89.9999);      /* Make sure the chart isn't being cast */
  127.   AA = MAX(AA, -89.9999);     /* on the precise north or south pole.  */
  128.   AA = DTOR(DecToDeg(AA));
  129.  
  130.   /* if parameter 'var' isn't set, then we can assume that the true time   */
  131.   /* has already been determined (as in a -rm switch time midpoint chart). */
  132.  
  133.   if (var) {
  134.     JD = (real)MdyToJulian(MM, DD, YY);
  135.     if (!progress || (operation & DASHp0) > 0)
  136.       T = ((JD-2415020.0)+TT/24.0-0.5)/36525.0;
  137.     else
  138.  
  139.       /* Determine actual time that a progressed chart is to be cast for. */
  140.  
  141.       T = (((Jdp-JD)/progday+JD)-2415020.0+TT/24.0-0.5)/36525.0;
  142.   }
  143.  
  144.   /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  145.   OB = DTOR(23.452294-0.0130125*T);
  146.  
  147.   Ln = Mod((933060-6962911*T+7.5*T*T)/3600.0);    /* Mean lunar node */
  148.   Off = (259205536.0*T+2013816.0)/3600.0;         /* Mean Sun        */
  149.   Off = 17.23*sin(DTOR(Ln))+1.27*sin(DTOR(Off))-(5025.64+1.11*T)*T;
  150.   Off = (Off-84038.27)/3600.0;
  151.   SD = ((operation & DASHs) ? Off : 0.0) + addfactor;
  152.   return Off;
  153. }
  154.  
  155.  
  156. /* Convert polar to rectangular coordinates. */
  157.  
  158. void PolToRec(A, R, X, Y)
  159. real A, R, *X, *Y;
  160. {
  161.   if (A == 0.0)
  162.     A = SMALL;
  163.   *X = R*cos(A);
  164.   *Y = R*sin(A);
  165. }
  166.  
  167.  
  168. /* Convert rectangular to polar coordinates. */
  169.  
  170. void RecToPol(X, Y, A, R)
  171. real X, Y, *A, *R;
  172. {
  173.   if (Y == 0.0)
  174.     Y = SMALL;
  175.   *R = sqrt(X*X+Y*Y);
  176.   *A = Angle(X, Y);
  177. }
  178.  
  179.  
  180. /* Convert rectangular to spherical coordinates. */
  181.  
  182. real RecToSph(B, L, O)
  183. real B, L, O;
  184. {
  185.   real R, Q, G;
  186.  
  187.   A = B; R = 1.0;
  188.   PolToRec(A, R, &X, &Y);
  189.   Q = Y; R = X; A = L;
  190.   PolToRec(A, R, &X, &Y);
  191.   G = X; X = Y; Y = Q;
  192.   RecToPol(X, Y, &A, &R);
  193.   A += O;
  194.   PolToRec(A, R, &X, &Y);
  195.   Q = ASIN(Y);
  196.   Y = X; X = G;
  197.   RecToPol(X, Y, &A, &R);
  198.   if (A < 0.0)
  199.     A += 2*PI;
  200.   G = A;
  201.   return G;
  202. }
  203.  
  204.  
  205. /* Do a coordinate transformation: Given a longitude and latitude value,    */
  206. /* return the new longitude and latitude values that the same location      */
  207. /* would have, were the equator tilted by a specified number of degrees.    */
  208. /* In other words, do a pole shift! This is used to convert among ecliptic, */
  209. /* equatorial, and local coordinates, each of which have zero declination   */
  210. /* in different planes. In other words, take into account the Earth's axis. */
  211.  
  212. void CoorXform(azi, alt, tilt)
  213. real *azi, *alt, tilt;
  214. {
  215.   real x, y, a1, l1;
  216.   real sinalt, cosalt, sinazi, sintilt, costilt;
  217.  
  218.   sinalt = sin(*alt); cosalt = cos(*alt); sinazi = sin(*azi);
  219.   sintilt = sin(tilt); costilt = cos(tilt);
  220.  
  221.   x = cosalt * sinazi * costilt;
  222.   y = sinalt * sintilt;
  223.   x -= y;
  224.   a1 = cosalt;
  225.   y = cosalt * cos(*azi);
  226.   l1 = Angle(y, x);
  227.   a1 = a1 * sinazi * sintilt + sinalt * costilt;
  228.   a1 = ASIN(a1);
  229.   *azi = l1; *alt = a1;
  230. }
  231.  
  232.  
  233. /* This is another subprocedure of CastChart(). Calculate a few variables */
  234. /* corresponding to the chart parameters that are used later on. The      */
  235. /* astrological vertex (object number twenty) is also calculated here.    */
  236.  
  237. void ComputeVariables()
  238. {
  239.   real R, R2, B, L, O, G;
  240.  
  241.   RA = DTOR(Mod((6.6460656+2400.0513*T+2.58E-5*T*T+TT)*15.0-OO));
  242.   R2 = RA; O = -OB; B = AA; A = R2; R = 1.0;
  243.   PolToRec(A, R, &X, &Y);
  244.   X *= cos(O);
  245.   RecToPol(X, Y, &A, &R);
  246.   MC = Mod(SD+RTOD(A));            /* Midheaven */
  247.   L = R2;
  248.   G = RecToSph(B, L, O);
  249. #if FALSE
  250.   Asc = Mod(SD+Mod(G+PI/2.0));     /* Ascendant */
  251. #endif
  252.   L= R2+PI; B = PI/2.0-dabs(B);
  253.   if (AA < 0.0)
  254.     B = -B;
  255.   G = RecToSph(B, L, O);
  256.   Vtx = Mod(SD+RTOD(G+PI/2.0));    /* Vertex */
  257. }
  258. #endif /* MATRIX */
  259.  
  260.  
  261. /*
  262. ******************************************************************************
  263. ** House Cusp Calculations.
  264. ******************************************************************************
  265. */
  266.  
  267. /* This is a subprocedure of HousePlace(). Given a zodiac position, return */
  268. /* which of the twelve houses it falls in. Remember that a special check   */
  269. /* has to be done for the house that spans 0 degrees Aries.                */
  270.  
  271. int HousePlaceIn(point)
  272. real point;
  273. {
  274.   int i = 0;
  275.  
  276.   point = Mod(point + 0.5/60.0/60.0);
  277.   do {
  278.     i++;
  279.   } while (!(i >= SIGNS ||
  280.       (point >= house[i] && point < house[Mod12(i+1)]) ||
  281.       (house[i] > house[Mod12(i+1)] &&
  282.       (point >= house[i] || point < house[Mod12(i+1)]))));
  283.   return i;
  284. }
  285.  
  286.  
  287. /* For each object in the chart, determine what house it belongs in. */
  288.  
  289. void HousePlace()
  290. {
  291.   int i;
  292.  
  293.   for (i = 1; i <= total; i++)
  294.     inhouse[i] = HousePlaceIn(planet[i]);
  295. }
  296.  
  297.  
  298. #ifdef MATRIX
  299. /* The following two functions calculate the midheaven and ascendant of  */
  300. /* the chart in question, based on time and location. They are also used */
  301. /* in some of the house cusp calculation routines as a quick way to get  */
  302. /* the 10th and 1st house cusps.                                         */
  303.  
  304. real CuspMidheaven()
  305. {
  306.   real MC;
  307.  
  308.   MC = ATAN(tan(RA)/cos(OB));
  309.   if (MC < 0.0)
  310.     MC += PI;
  311.   if (RA > PI)
  312.     MC += PI;
  313.   return Mod(RTOD(MC)+SD);
  314. }
  315.  
  316. real CuspAscendant()
  317. {
  318.   real Asc;
  319.  
  320.   Asc = Angle(-sin(RA)*cos(OB)-tan(AA)*sin(OB), cos(RA));
  321.   return Mod(RTOD(Asc)+SD);
  322. }
  323.  
  324.  
  325. /* These are various different algorithms for calculating the house cusps: */
  326.  
  327. real CuspPlacidus(deg, FF, neg)
  328. real deg, FF;
  329. bool neg;
  330. {
  331.   real LO, R1, R2, XS;
  332.   int i;
  333.  
  334.   R1 = RA+DTOR(deg);
  335.   if (neg)
  336.     X = 1.0;
  337.   else
  338.     X = -1.0;
  339.   for (i = 1; i <= 10; i++) {
  340.  
  341.     /* This formula works except at 0 latitude (AA == 0.0). */
  342.  
  343.     XS = X*sin(R1)*tan(OB)*tan(AA == 0.0 ? 0.0001 : AA);
  344.     XS = ACOS(XS);
  345.     if (XS < 0.0)
  346.       XS += PI;
  347.     if (neg)
  348.       R2 = RA+PI-(XS/FF);
  349.     else
  350.       R2 = RA+(XS/FF);
  351.     R1 = R2;
  352.   }
  353.   LO = ATAN(tan(R1)/cos(OB));
  354.   if (LO < 0.0)
  355.     LO += PI;
  356.   if (sin(R1) < 0.0)
  357.     LO += PI;
  358.   return RTOD(LO);
  359. }
  360.  
  361. void HousePlacidus()
  362. {
  363.   int i;
  364.  
  365.   house[1] = Mod(Asc-SD);
  366.   house[4] = Mod(MC+DEGHALF-SD);
  367.   house[5] = CuspPlacidus(30.0, 3.0, FALSE) + DEGHALF;
  368.   house[6] = CuspPlacidus(60.0, 1.5, FALSE) + DEGHALF;
  369.   house[2] = CuspPlacidus(120.0, 1.5, TRUE);
  370.   house[3] = CuspPlacidus(150.0, 3.0, TRUE);
  371.   for (i = 1; i <= SIGNS; i++) {
  372.     if (i <= 6)
  373.       house[i] = Mod(house[i]+SD);
  374.     else
  375.       house[i] = Mod(house[i-6]+DEGHALF);
  376.   }
  377. }
  378.  
  379. void HouseKoch()
  380. {
  381.   real A1, A2, A3, KN, D;
  382.   int i;
  383.  
  384.   A1 = sin(RA)*tan(AA)*tan(OB);
  385.   A1 = ASIN(A1);
  386.   for (i = 1; i <= SIGNS; i++) {
  387.     D = Mod(60.0+30.0*(real)i);
  388.     A2 = D/DEGQUAD-1.0; KN = 1.0;
  389.     if (D >= DEGHALF) {
  390.       KN = -1.0;
  391.       A2 = D/DEGQUAD-3.0;
  392.     }
  393.     A3 = DTOR(Mod(RTOD(RA)+D+A2*RTOD(A1)));
  394.     X = Angle(cos(A3)*cos(OB)-KN*tan(AA)*sin(OB), sin(A3));
  395.     house[i] = Mod(RTOD(X)+SD);
  396.   }
  397. }
  398.  
  399. void HouseEqual()
  400. {
  401.   int i;
  402.  
  403.   for (i = 1; i <= SIGNS; i++)
  404.     house[i] = Mod(Asc-30.0+30.0*(real)i);
  405. }
  406.  
  407. void HouseCampanus()
  408. {
  409.   real KO, DN;
  410.   int i;
  411.  
  412.   for (i = 1; i <= SIGNS; i++) {
  413.     KO = DTOR(60.000001+30.0*(real)i);
  414.     DN = ATAN(tan(KO)*cos(AA));
  415.     if (DN < 0.0)
  416.       DN += PI;
  417.     if (sin(KO) < 0.0)
  418.       DN += PI;
  419.     X = Angle(cos(RA+DN)*cos(OB)-sin(DN)*tan(AA)*sin(OB), sin(RA+DN));
  420.     house[i] = Mod(RTOD(X)+SD);
  421.   }
  422. }
  423.  
  424. void HouseMeridian()
  425. {
  426.   real D;
  427.   int i;
  428.  
  429.   for (i = 1; i <= SIGNS; i++) {
  430.     D = DTOR(60.0+30.0*(real)i);
  431.     X = Angle(cos(RA+D)*cos(OB), sin(RA+D));
  432.     house[i] = Mod(RTOD(X)+SD);
  433.   }
  434. }
  435.  
  436. void HouseRegiomontanus()
  437. {
  438.   real D;
  439.   int i;
  440.  
  441.   for (i = 1; i <= SIGNS; i++) {
  442.     D = DTOR(60.0+30.0*(real)i);
  443.     X = Angle(cos(RA+D)*cos(OB)-sin(D)*tan(AA)*sin(OB), sin(RA+D));
  444.     house[i] = Mod(RTOD(X)+SD);
  445.   }
  446. }
  447.  
  448. void HousePorphyry()
  449. {
  450.   int i;
  451.  
  452.   X = Asc-MC;
  453.   if (X < 0.0)
  454.     X += DEGREES;
  455.   Y = X/3.0;
  456.   for (i = 1; i <= 2; i++)
  457.     house[i+4] = Mod(DEGHALF+MC+i*Y);
  458.   X = Mod(DEGHALF+MC)-Asc;
  459.   if (X < 0.0)
  460.     X += DEGREES;
  461.   house[1]=Asc;
  462.   Y = X/3.0;
  463.   for (i = 1; i <= 3; i++)
  464.     house[i+1] = Mod(Asc+i*Y);
  465.   for (i = 1; i <= 6; i++)
  466.     house[i+6] = Mod(house[i]+DEGHALF);
  467. }
  468.  
  469. void HouseMorinus()
  470. {
  471.   real D;
  472.   int i;
  473.  
  474.   for (i = 1; i <= SIGNS; i++) {
  475.     D = DTOR(60.0+30.0*(real)i);
  476.     X = Angle(cos(RA+D), sin(RA+D)*cos(OB));
  477.     house[i] = Mod(RTOD(X)+SD);
  478.   }
  479. }
  480.  
  481. real CuspTopocentric(deg)
  482. real deg;
  483. {
  484.   real OA, X, LO;
  485.  
  486.   OA = Mod(RA+DTOR(deg));
  487.   X = ATAN(tan(AA)/cos(OA));
  488.   LO = ATAN(cos(X)*tan(OA)/cos(X+OB));
  489.   if (LO < 0.0)
  490.     LO += PI;
  491.   if (sin(OA) < 0.0)
  492.     LO += PI;
  493.   return LO;
  494. }
  495.  
  496. void HouseTopocentric()
  497. {
  498.   real TL, P1, P2, LT;
  499.   int i;
  500.  
  501.   modulus = 2.0*PI;
  502.   house[4] = Mod(DTOR(MC+DEGHALF-SD));
  503.   TL = tan(AA); P1 = ATAN(TL/3.0); P2 = ATAN(TL/1.5); LT = AA;
  504.   AA = P1; house[5] = CuspTopocentric(30.0) + PI;
  505.   AA = P2; house[6] = CuspTopocentric(60.0) + PI;
  506.   AA = LT; house[1] = CuspTopocentric(90.0);
  507.   AA = P2; house[2] = CuspTopocentric(120.0);
  508.   AA = P1; house[3] = CuspTopocentric(150.0);
  509.   AA = LT; modulus = DEGREES;
  510.   for (i = 1; i <= 6; i++) {
  511.     house[i] = Mod(RTOD(house[i])+SD);
  512.     house[i+6] = Mod(house[i]+DEGHALF);
  513.   }
  514. }
  515. #endif /* MATRIX */
  516.  
  517. /* This house system is just like the Equal system except that we start */
  518. /* our 12 equal segments from the Midheaven instead of the Ascendant.   */
  519.  
  520. void HouseEqualMidheaven()
  521. {
  522.   int i;
  523.  
  524.   for (i = 1; i <= SIGNS; i++)
  525.     house[i] = Mod(MC-270.0+30.0*(real)(i-1));
  526. }
  527.  
  528. /* This is a new house system similar in philosophy to Porphyry houses.   */
  529. /* Instead of just trisecting the difference in each quadrant, we do a    */
  530. /* smooth sinusoidal distribution of the difference around all the cusps. */
  531.  
  532. void HousePorphyryNeo()
  533. {
  534.   real delta;
  535.   int i;
  536.  
  537.   delta = (MinDistance(MC, Asc) - DEGQUAD)/4.0;
  538.   house[_LIB] = Mod(Asc+DEGHALF); house[_CAP] = MC;
  539.   house[_AQU] = Mod(house[_CAP] + 30.0 + delta   + SD);
  540.   house[_PIS] = Mod(house[_AQU] + 30.0 + delta*2 + SD);
  541.   house[_SAG] = Mod(house[_CAP] - 30.0 + delta   + SD);
  542.   house[_SCO] = Mod(house[_SAG] - 30.0 + delta*2 + SD);
  543.   for (i = _ARI; i < _LIB; i++)
  544.     house[i] = Mod(house[i+6]-DEGHALF);
  545. }
  546.  
  547. /* In "null" houses, the cusps are always fixed to start at their            */
  548. /* corresponding sign, i.e. the 1st house is always at 0 degrees Aries, etc. */
  549.  
  550. void HouseNull()
  551. {
  552.   int i;
  553.  
  554.   for (i = 1; i <= SIGNS; i++)
  555.     house[i] = Mod(STOZ(i)+SD);
  556. }
  557.  
  558.  
  559. /* Calculate the house cusp positions, using the specified algorithm. */
  560.  
  561. void ComputeHouses(housesystem)
  562. int housesystem;
  563. {
  564.   char string[STRING];
  565.  
  566.   if (dabs(AA) > DTOR(DEGQUAD-TROPIC) && housesystem < 2) {
  567.     sprintf(string,
  568.       "The %s system of houses is not defined at extreme latitudes.",
  569.       systemname[housesystem]);
  570.     PrintError(string);
  571.     Terminate(_FATAL);
  572.   }
  573.   switch (housesystem) {
  574.   case  1: HouseKoch();           break;
  575.   case  2: HouseEqual();          break;
  576.   case  3: HouseCampanus();       break;
  577.   case  4: HouseMeridian();       break;
  578.   case  5: HouseRegiomontanus();  break;
  579.   case  6: HousePorphyry();       break;
  580.   case  7: HouseMorinus();        break;
  581.   case  8: HouseTopocentric();    break;
  582.   case  9: HouseEqualMidheaven(); break;
  583.   case 10: HousePorphyryNeo();    break;
  584.   case 11: HouseNull();           break;
  585.   default: HousePlacidus();
  586.   }
  587. }
  588.  
  589.  
  590. #ifdef MATRIX
  591. /*
  592. ******************************************************************************
  593. ** Planetary Position Calculations.
  594. ******************************************************************************
  595. */
  596.  
  597. /* Read the next three values from the planet data stream, and return them */
  598. /* combined as the coefficients of a quadratic equation in the chart time. */
  599.  
  600. real ReadThree()
  601. {
  602.   real S0, S1, S2;
  603.  
  604.   S0 = ReadPlanetData(); S1 = ReadPlanetData();
  605.   S2 = ReadPlanetData();
  606.   return S0 = DTOR(S0+S1*T+S2*T*T);
  607. }
  608.  
  609.  
  610. /* Another coordinate transformation. This is used by the ComputePlanets() */
  611. /* procedure to rotate rectangular coordinates by a certain amount.        */
  612.  
  613. real RecToSph2(AP, AN, IN)
  614. real AP, AN, IN;
  615. {
  616.   real R, D, G;
  617.  
  618.   RecToPol(X, Y, &A, &R); A += AP; PolToRec(A, R, &X, &Y);
  619.   D = X; X = Y; Y = 0.0; RecToPol(X, Y, &A, &R);
  620.   A += IN; PolToRec(A, R, &X, &Y);
  621.   G = Y; Y = X; X = D; RecToPol(X, Y, &A, &R); A += AN;
  622.   if (A < 0.0)
  623.     A += 2.0*PI;
  624.   PolToRec(A, R, &X, &Y);
  625.   return G;
  626. }
  627.  
  628.  
  629. /* Calculate some harmonic delta error correction factors to add onto the */
  630. /* coordinates of Jupiter through Pluto, for better accuracy.             */
  631.  
  632. void ErrorCorrect(ind, x, y, z)
  633. int ind;
  634. real *x, *y, *z;
  635. {
  636.   real U, V, W, S0, T0[4];
  637.   int IK, IJ, errorindex;
  638.  
  639.   errorindex = errorcount[ind];
  640.   for (IK = 1; IK <= 3; IK++) {
  641.     if (ind == 6 && IK == 3) {
  642.       T0[3] = 0.0;
  643.       break;
  644.     }
  645.     if (IK == 3)
  646.       errorindex--;
  647.     S0 = ReadThree(); A = 0.0;
  648.     for (IJ = 1; IJ <= errorindex; IJ++) {
  649.       U = ReadPlanetData(); V = ReadPlanetData(); W = ReadPlanetData();
  650.       A = A+DTOR(U)*cos((V*T+W)*PI/DEGHALF);
  651.     }
  652.     T0[IK] = RTOD(S0+A);
  653.   }
  654.   *x += T0[2]; *y += T0[1]; *z += T0[3];
  655. }
  656.  
  657.  
  658. /* Another subprocedure of the ComputePlanets() routine. Convert the final */
  659. /* rectangular coordinates of a planet to zodiac position and declination. */
  660.  
  661. void ProcessPlanet(ind, aber)
  662. int ind;
  663. real aber;
  664. {
  665.   real ang, rad;
  666.  
  667.   RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  668.   planet[ind] = Mod(RTOD(ang) /*+ NU*/ - aber + SD);
  669.   RecToPol(rad, spacez[ind], &ang, &rad);
  670.   if (centerplanet == _SUN && ind == _SUN)
  671.     ang = 0.0;
  672.   planetalt[ind] = RTOD(ang);
  673. }
  674.  
  675.  
  676. /* This is probably the heart of the whole program of Astrolog. Calculate  */
  677. /* the position of each body that orbits the Sun. A heliocentric chart is  */
  678. /* most natural; extra calculation is needed to have other central bodies. */
  679.  
  680. void ComputePlanets()
  681. {
  682.   real helio[BASE+1], helioalt[BASE+1], helioret[BASE+1],
  683.     heliox[BASE+1], helioy[BASE+1], helioz[BASE+1];
  684.   real aber = 0.0, AU, E, EA, E1, M, XW, YW, AP, AN, IN, G, XS, YS, ZS;
  685.   int ind = 1, i;
  686.  
  687.   datapointer = planetdata;
  688.   while (ind <= (operation & DASHu ? BASE : PLANETS+1)) {
  689.     modulus = 2.0*PI;
  690.     EA = M = Mod(ReadThree());      /* Calculate mean anomaly */
  691.     E = RTOD(ReadThree());          /* Calculate eccentricity */
  692.     for (i = 1; i <= 5; i++)
  693.       EA = M+E*sin(EA);             /* Solve Kepler's equation */
  694.     AU = ReadPlanetData();          /* Semi-major axis         */
  695.     E1 = 0.01720209/(pow(AU, 1.5)*
  696.       (1.0-E*cos(EA)));                     /* Begin velocity coordinates */
  697.     XW = -AU*E1*sin(EA);                    /* Perifocal coordinates      */
  698.     YW = AU*E1*pow(1.0-E*E,0.5)*cos(EA);
  699.     AP = ReadThree(); AN = ReadThree();
  700.     IN = ReadThree();                           /* Calculate inclination  */
  701.     X = XW; Y = YW; G = RecToSph2(AP, AN, IN);  /* Rotate velocity coords */
  702.     heliox[ind] = X; helioy[ind] = Y;
  703.     helioz[ind] = G;                       /* Helio ecliptic rectangtular */
  704.     modulus = DEGREES;
  705.     X = AU*(cos(EA)-E);                 /* Perifocal coordinates for        */
  706.     Y = AU*sin(EA)*pow(1.0-E*E,0.5);    /* rectangular position coordinates */
  707.     G = RecToSph2(AP, AN, IN);          /* Rotate for rectangular */
  708.     XS = X; YS = Y; ZS = G;             /* position coordinates   */
  709.     if (ind >= _JUP && ind <= _PLU)
  710.       ErrorCorrect(ind, &XS, &YS, &ZS);
  711.     ret[ind] =                                        /* Helio daily motion */
  712.       (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
  713.     spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
  714.     ProcessPlanet(ind, 0.0);
  715.     ind += (ind == 1 ? 2 : (ind != PLANETS+1 ? 1 : 10));
  716.   }
  717.   spacex[0] = spacey[0] = spacez[0] = 0.0;
  718.  
  719.   if (!centerplanet) {
  720.     if (exdisplay & DASHv0)
  721.       for (i = 0; i <= BASE; i++)    /* Use relative velocity */
  722.         ret[i] = DTOR(1.0);          /* if -v0 is in effect   */
  723.     return;
  724.   }
  725. #endif /* MATRIX */
  726.  
  727.   /* A second loop is needed for geocentric charts or central bodies other */
  728.   /* than the Sun. For example, we can't find the position of Mercury in   */
  729.   /* relation to Pluto until we know the position of Pluto in relation to  */
  730.   /* the Sun, and since Mercury is calculated first, another pass needed.  */
  731.  
  732.   ind = centerplanet;
  733.   for (i = 0; i <= BASE; i++) {
  734.     helio[i]    = planet[i];
  735.     helioalt[i] = planetalt[i];
  736.     helioret[i] = ret[i];
  737.     if (i != _MOO && i != ind) {
  738.       spacex[i] -= spacex[ind];
  739.       spacey[i] -= spacey[ind];
  740.       spacez[i] -= spacez[ind];
  741.     }
  742.   }
  743.   spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  744.   SwapReal(&spacex[0], &spacex[ind]);
  745.   SwapReal(&spacey[0], &spacey[ind]);    /* Do some swapping - we want   */
  746.   SwapReal(&spacez[0], &spacez[ind]);    /* the central body to be in    */
  747.   SwapReal(&spacex[1], &spacex[ind]);    /* object position number zero. */
  748.   SwapReal(&spacey[1], &spacey[ind]);
  749.   SwapReal(&spacez[1], &spacez[ind]);
  750.   for (i = 1; i <= (operation & DASHu ? BASE : PLANETS+1);
  751.       i += (i == 1 ? 2 : (i != PLANETS+1 ? 1 : 10))) {
  752.     XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
  753.     if (ind != _SUN || i != _SUN)
  754.       ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
  755.         (XS*XS+YS*YS);
  756.     if (ind == _SUN)
  757.       aber = 0.0057756*sqrt(XS*XS+YS*YS+ZS*ZS)*RTOD(ret[i]);  /* Aberration */
  758.     ProcessPlanet(i, aber);
  759.     if (exdisplay & DASHv0)                 /* Use relative velocity */
  760.       ret[i] = DTOR(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  761.   }
  762. }
  763.  
  764.  
  765. #ifdef MATRIX
  766. /*
  767. ******************************************************************************
  768. ** Lunar Position Calculations
  769. ******************************************************************************
  770. */
  771.  
  772. /* Calculate the position and declination of the Moon, and the Moon's North  */
  773. /* Node. This has to be done separately from the other planets, because they */
  774. /* all orbit the Sun, while the Moon orbits the Earth.                       */
  775.  
  776. void ComputeLunar(moonlo, moonla, nodelo, nodela)
  777. real *moonlo, *moonla, *nodelo, *nodela;
  778. {
  779.   real LL, G, N, G1, D, L, ML, L1, MB, T1, M = 3600.0;
  780.  
  781.   LL = 973563.0+1732564379.0*T-4.0*T*T;  /* Compute mean lunar longitude     */
  782.   G = 1012395.0+6189.0*T;                /* Sun's mean longitude of perigee  */
  783.   N = 933060.0-6962911.0*T+7.5*T*T;      /* Compute mean lunar node          */
  784.   G1 = 1203586.0+14648523.0*T-37.0*T*T;  /* Mean longitude of lunar perigee  */
  785.   D = 1262655.0+1602961611.0*T-5.0*T*T;  /* Mean elongation of Moon from Sun */
  786.   L = (LL-G1)/M; L1 = ((LL-D)-G)/M;      /* Some auxiliary angles            */
  787.   T1 = (LL-N)/M; D = D/M; Y = 2.0*D;
  788.  
  789.   /* Compute Moon's perturbations. */
  790.  
  791.   ML = 22639.6*SIND(L)-4586.4*SIND(L-Y)+2369.9*SIND(Y)+769.0*SIND(2.0*L)-
  792.     669.0*SIND(L1)-411.6*SIND(2.0*T1)-212.0*SIND(2.0*L-Y)-206.0*SIND(L+L1-Y);
  793.   ML += 192.0*SIND(L+Y)-165.0*SIND(L1-Y)+148.0*SIND(L-L1)-125.0*SIND(D)-110.0*
  794.     SIND(L+L1)-55.0*SIND(2.0*T1-Y)-45.0*SIND(L+2.0*T1)+40.0*SIND(L-2.0*T1);
  795.  
  796.   *moonlo = G = Mod((LL+ML)/M+SD);       /* Lunar longitude */
  797.  
  798.   /* Compute lunar latitude. */
  799.  
  800.   MB = 18461.5*SIND(T1)+1010.0*SIND(L+T1)-999.0*SIND(T1-L)-624.0*SIND(T1-Y)+
  801.     199.0*SIND(T1+Y-L)-167.0*SIND(L+T1-Y);
  802.   MB += 117.0*SIND(T1+Y)+62.0*SIND(2.0*L+T1)-
  803.     33.0*SIND(T1-Y-L)-32.0*SIND(T1-2.0*L)-30.0*SIND(L1+T1-Y);
  804.   *moonla = MB =
  805.     Sgn(MB)*((dabs(MB)/M)/DEGREES-floor((dabs(MB)/M)/DEGREES))*DEGREES;
  806.  
  807.   /* Compute position of the north node. One can compute the true North */
  808.   /* Node here if they like, but the Mean Node seems to be the one used */
  809.   /* in astrology, so that's the one Astrolog returns by default.       */
  810.  
  811. #ifdef TRUENODE
  812.   N = N+5392.0*SIND(2.0*T1-Y)-541.0*SIND(L1)-442.0*SIND(Y)+423.0*SIND(2.0*T1)-
  813.     291.0*SIND(2.0*L-2.0*T1);
  814. #endif
  815.   *nodelo = Mod(N/M+SD);
  816.   *nodela = 0.0;
  817. }
  818. #endif /* MATRIX */
  819.  
  820.  
  821. /*
  822. ******************************************************************************
  823. ** Star Position Calculations
  824. ******************************************************************************
  825. */
  826.  
  827. /* This is used by the chart calculation routine to calculate the positions */
  828. /* of the fixed stars. Since the stars don't move in the sky over time,     */
  829. /* getting their positions is mostly just reading info from an array and    */
  830. /* converting it to the correct reference frame. However, we have to add    */
  831. /* in the correct precession for the tropical zodiac, and sort the final    */
  832. /* index list based on what order the stars are supposed to be printed in.  */
  833.  
  834. void ComputeStars(SD)
  835. real SD;
  836. {
  837.   int i, j;
  838.   real x, y, z;
  839.  
  840.   /* Read in star positions. */
  841.  
  842.   for (i = 1; i <= STARS; i++) {
  843.     x = stardata[i*6-6]; y = stardata[i*6-5]; z = stardata[i*6-4];
  844.     planet[BASE+i] = DTOR(x*DEGREES/24.0+y*15.0/60.0+z*0.25/60.0);
  845.     x = stardata[i*6-3]; y = stardata[i*6-2]; z = stardata[i*6-1];
  846.     planetalt[BASE+i] = DTOR(x+y/60.0+z/60.0/60.0);
  847.     EquToEcl(&planet[BASE+i], &planetalt[BASE+i]);           /* Convert to */
  848.     planet[BASE+i] = Mod(RTOD(planet[BASE+i])+SD2000+SD);    /* ecliptic.  */
  849.     planetalt[BASE+i] = RTOD(planetalt[BASE+i]);
  850.     starname[i] = i;
  851.   }
  852.  
  853.   /* Sort the index list if -Uz, -Ul, -Un, or -Ub switch in effect. */
  854.  
  855.   if (universe > 1) for (i = 2; i <= STARS; i++) {
  856.     j = i-1;
  857.  
  858.     /* Compare star names for -Un switch. */
  859.  
  860.     if (universe == 'n') while (j > 0 && StringCmp(
  861.       objectname[BASE+starname[j]], objectname[BASE+starname[j+1]]) > 0) {
  862.       SWAP(starname[j], starname[j+1]);
  863.       j--;
  864.  
  865.     /* Compare star brightnesses for -Ub switch. */
  866.  
  867.     } else if (universe == 'b') while (j > 0 &&
  868.       starbright[starname[j]] > starbright[starname[j+1]]) {
  869.       SWAP(starname[j], starname[j+1]);
  870.       j--;
  871.  
  872.     /* Compare star zodiac locations for -Uz switch. */
  873.  
  874.     } else if (universe == 'z') while (j > 0 &&
  875.       planet[BASE+starname[j]] > planet[BASE+starname[j+1]]) {
  876.       SWAP(starname[j], starname[j+1]);
  877.       j--;
  878.  
  879.     /* Compare star declinations for -Ul switch. */
  880.  
  881.     } else if (universe == 'l') while (j > 0 &&
  882.       planetalt[BASE+starname[j]] < planetalt[BASE+starname[j+1]]) {
  883.       SWAP(starname[j], starname[j+1]);
  884.       j--;
  885.     }
  886.   }
  887. }
  888.  
  889.  
  890. /*
  891. ******************************************************************************
  892. ** Chart Calculation.
  893. ******************************************************************************
  894. */
  895.  
  896. #ifdef PLACALC
  897. /* Compute the positions of the planets at a certain time using the Placalc */
  898. /* accurate formulas and ephemeris. This will superseed the Matrix routine  */
  899. /* values and is only called with the -b switch is in effect. Not all       */
  900. /* objects or modes are available using this, but some additional values    */
  901. /* such as Moon and Node velocities not available without -b are. (This is  */
  902. /* the one place in Astrolog which calls the Placalc package functions.)    */
  903.  
  904. void ComputePlacalc(t)
  905. real t;
  906. {
  907.   int i;
  908.   real r1, r2, r3, r4, r;
  909.  
  910.   /* We can compute the positions of Sun through Pluto, Chiron, and the */
  911.   /* North Node using Placalc. The other object must be done elsewhere. */
  912.  
  913.   for (i = _SUN; i <= _NOD; i++) if (i <= _CHI || i >= _NOD) {
  914.     if (PlacalcPlanet(i, t*36525.0+2415020.0, centerplanet != _SUN,
  915.       &r1, &r2, &r3, &r4)) {
  916.  
  917.       /* Note that this can't compute charts with central planets other */
  918.       /* than the Sun or Earth or relative velocities in current state. */
  919.  
  920.       planet[i]    = Mod(r1 + SD);
  921.       planetalt[i] = r2;
  922.       ret[i]       = DTOR(r3);
  923.  
  924.       /* Compute x,y,z coordinates from azimuth, altitude, and distance. */
  925.  
  926.       spacez[i] = r4*SIND(planetalt[i]);
  927.       r         = r4*COSD(planetalt[i]);
  928.       spacex[i] = r*COSD(planet[i]);
  929.       spacey[i] = r*SIND(planet[i]);
  930.     }
  931.   }
  932. }
  933. #endif
  934.  
  935.  
  936. /* This is probably the main routine in all of Astrolog. It generates a   */
  937. /* chart, calculating the positions of all the celestial bodies and house */
  938. /* cusps, based on the current chart information, and saves them for use  */
  939. /* by any of the display routines.                                        */
  940.  
  941. real CastChart(var)
  942. int var;
  943. {
  944.   int i, k;
  945.   real housetemp[SIGNS+1], Off = 0.0, j;
  946.  
  947.   if (MM == -1) {
  948.  
  949.     /* Hack: If month is negative, then we know chart was read in through a  */
  950.     /* -o0 position file, so the planet positions are already in the arrays. */
  951.  
  952.     MC = planet[_MC]; Asc = planet[_ASC]; Vtx = planet[_VTX];
  953.   } else {
  954.     Off = ProcessInput(var);
  955.     ComputeVariables();
  956.     if (operation & DASHG)          /* Check for -G geodetic chart. */
  957.       RA = DTOR(Mod(-OO));
  958.     MC  = CuspMidheaven();          /* Calculate our Ascendant & Midheaven. */
  959.     Asc = CuspAscendant();
  960.     ComputeHouses(housesystem);     /* Go calculate house cusps. */
  961.     for (i = 1; i <= total; i++) {
  962.       planetalt[i] = 0.0;           /* Assume on ecliptic unless we say so.  */
  963.       ret[i] = 1.0;                 /* Assume direct until we say otherwise. */
  964.     }
  965.  
  966.     /* Go calculate planet, Moon, and North Node positions. */
  967.  
  968.     ComputePlanets();
  969.     ComputeLunar(&planet[_MOO], &planetalt[_MOO],
  970.       &planet[_NOD], &planetalt[_NOD]);
  971.     ret[_NOD] = -1.0;
  972.  
  973.     /* Compute more accurate ephemeris positions for certain objects. */
  974.  
  975. #ifdef PLACALC
  976.     if (placalc)
  977.       ComputePlacalc(T);
  978. #endif
  979.  
  980.     /* Calculate position of Part of Fortune. */
  981.  
  982.     j = planet[_MOO]-planet[_SUN];
  983.     j = dabs(j) < DEGQUAD ? j : j - Sgn(j)*DEGREES;
  984.     planet[_FOR] = Mod(j+Asc);
  985.  
  986.     /* Fill in "planet" positions corresponding to house cusps. */
  987.  
  988.     planet[_MC] = MC; planet[_ASC] = Asc; planet[_VTX] = Vtx;
  989.     planet[C_LO]   = house[11]; planet[C_LO+1] = house[12];
  990.     planet[C_LO+2] = house[2];  planet[C_LO+3] = house[3];
  991.   }
  992.  
  993.   /* Go calculate star positions if -U switch in effect. */
  994.  
  995.   if (universe)
  996.     ComputeStars(operation & DASHs ? 0.0 : -Off);
  997.  
  998.   /* Now, we may have to modify the base positions we calculated above based */
  999.   /* on what type of chart we are generating.                                */
  1000.  
  1001.   if (operation & DASHp0) {         /* Are we doing a -p0 solar arc chart?   */
  1002.     for (i = 1; i <= total; i++)
  1003.       planet[i] = Mod(planet[i] + (Jdp - JD) / progday);
  1004.     for (i = 1; i <= SIGNS; i++)
  1005.       house[i]  = Mod(house[i]  + (Jdp - JD) / progday);
  1006.     }
  1007.   if (multiplyfactor > 1)           /* Are we doing a -x harmonic chart?     */
  1008.     for (i = 1; i <= total; i++)
  1009.       planet[i] = Mod(planet[i] * (real)multiplyfactor);
  1010.   if (onasc) {
  1011.     if (onasc > 0)                  /* Is -1 put on Ascendant in effect?     */
  1012.       j = planet[onasc]-Asc;
  1013.     else                            /* Or -2 put object on Midheaven switch? */
  1014.       j = planet[-onasc]-MC;
  1015.     for (i = 1; i <= SIGNS; i++)    /* If so, rotate the houses accordingly. */
  1016.       house[i] = Mod(house[i]+j);
  1017.   }
  1018.  
  1019.   /* Check to see if we are -F forcing any objects to be particular values. */
  1020.  
  1021.   for (i = 1; i <= total; i++)
  1022.     if (force[i] != 0.0) {
  1023.       planet[i] = force[i]-DEGREES;
  1024.       planetalt[i] = ret[i] = 0.0;
  1025.     }
  1026.   HousePlace();               /* Figure out what house everything falls in. */
  1027.  
  1028.   /* If -f domal chart switch in effect, switch planet and house positions. */
  1029.  
  1030.   if (operation & DASHf) {
  1031.     for (i = 1; i <= total; i++) {
  1032.       k = inhouse[i];
  1033.       inhouse[i] = ZTOS(planet[i]);
  1034.       planet[i] = STOZ(k)+MinDistance(house[k], planet[i]) /
  1035.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  1036.     }
  1037.     for (i = 1; i <= SIGNS; i++) {
  1038.       k = HousePlaceIn(STOZ(i));
  1039.       housetemp[i] = STOZ(k)+MinDistance(house[k], STOZ(i)) /
  1040.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  1041.     }
  1042.     for (i = 1; i <= SIGNS; i++)
  1043.       house[i] = housetemp[i];
  1044.   }
  1045.  
  1046.   /* If -3 decan chart switch in effect, edit planet positions accordingly. */
  1047.  
  1048.   if (operation & DASH3)
  1049.     for (i = 1; i <= total; i++) {
  1050.       k = ZTOS(planet[i]);
  1051.       j = planet[i] - STOZ(k);
  1052.       k = Mod12(k + 4*((int)floor(j/10.0)));
  1053.       j = (j - floor(j/10.0)*10.0)*3.0;
  1054.       planet[i] = STOZ(k)+j;
  1055.       HousePlace();
  1056.     }
  1057.   return T;
  1058. }
  1059.  
  1060. /* formulas.c */
  1061.